STATISTICAL  EXTRAPOLATION 


Dennis  Randolph  Oldson 


United    States 
Naval   Postgraduate  School 


T 


STATISTICAL  EXTRAPOLATION 


by 


Dennis  Randolph  Oldson 


Thesis  Advisor: 


Sydney  R.  Parker 


September  1971 


Approved  {±o\  puhtlc  n.a  lease.;   cUAtiibution  Luitlmitzd. 


Statistical  Extrapolation 


by 


Dennis  Randolph  Oldson 

Lieutenant,  United  States  Navy 

B.S.E.E.,  University  of  Colorado,  1965 


Submitted  in  partial  fulfillment  of  the 
requirements  for  the  degree  of 


ELECTRICAL  ENGINEER 


from  the 

NAVAL  POSTGRADUATE  SCHOOL 
September  1971 


ABSTRACT 

The  mathematics  of  extrapolating  known  statistics  of 
components  to  the  probability  density  function  of  a  system's 
performance  measure  is  considered.   Quadrature  sum  inte- 
gration schemes  for  evaluating  the  resulting  required  inte- 
gration are  examined,  and  alternate  integral  approximation 
schemes  are  developed  utilizing  Monte  Carlo  methods.   A 
simple  electrical  circuit  example  illustrates  the  use  of 
these  techniques. 
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I.   INTRODUCTION 

The  desirability  of  being  able  to  statistically  predict 
a  system's  performance  measure  is  obvious  when  one  considers 
the  mass  production  of  systems.   Since  performance  measures 
of  certain  values  may  not  be  desirable,  it  is  generally  un- 
profitable to  produce  such  systems  when  the  probability  that 
these  undesirable  values  can  exist  is  beyond  the  limits 
determined  by  the  requirements  of  the  situation. 

Presently  performance  data  for  systems  are  obtained  in 
several  ways:   (1)   A  worst  case  analysis  [1]  to  determine  if 
undesirable  performance  will  occur  within  the  parameter  toler- 
ance; (2)   A  moment  method  [1]  in  which  the  performance 
measure  is  sampled  so  that  its  mean,  variance,  and  higher 
order  moments  are  determined  for  approximation  by  known  dis- 
tribution functions;  (3)   A  functional  formulas  method  [2]  in 
which  performance  measure's  distributions  are  obtained  from 
breaking  performance  measures  into  a  series  of  products  and 
sums  for  individual  integration;  (4)   Monte  Carlo  methods 
[3,  10]  in  which  performance  measure  distributions  are  gener- 
ated  by  creation  of  random  processes  by  which  the  parameters 
of  the  random  process  lead  to  an  approximation  to  the  desired 
distribution. 

The  purpose  of  this  thesis  is  to  investigate  the  mathe- 
matics required  in  obtaining  a  performance  measure  and  at- 
tempt to  arrive  at  a  reasonably  fast  and  accurate  method  of 


extrapolating  the  statistics  of  a  system's  component  variables 
to  create  the  distribution  of  the  desired  performance  measure. 

Chapter  II  examines  the  mathematical  form  and  statistical 
nature  of  the  performance  measure  for  implicit  and  explicit 
functions  and  presents  two  digital  computer  integration 
schemes  using  Gauss  quadrature  sums  for  evaluating  the  re- 
sulting required  integration. 

Chapter  III  investigates  crude  Monte  Carlo  methods  for 
determining  performance  measures,  and  the  confidence  limits 
involved  in  the  method.   Then  Monte  Carlo  schemes  to  estimate 
single  integrals  are  examined.   The  multiple  integral  exten- 
sions of  the  single  integral  estimates  are  then  derived. 

Chapter  IV  discusses  the  relative  merits  of  the  methods 
outlined  in  Chapters  II  and  III,  and  gives  a  comparison  of 
the  results  of  these  methods  applied  to  a  simple  example  con- 
sisting of  an  electrical  circuit  in  which  the  component 
values  are  the  random  variables. 


II.   STATISTICAL  ANALYSIS 

The  probability  density  function  of  a  system's  perform- 
ance measure  can  be  determined  mathematically  from  knowledge 
of  the  functional  relation  between  the  performance  measure 
and  the  system  component-variables,  and  knowledge  of  the 
joint  probability  density  function  of  the  component- 
variables.   Most  often  this  requires  integrating  a  complex 
function  over  several  variables.   Although  the  integration 
can  be  performed  analytically  in  some  cases,  digital  computer 
integration  schemes  must  normally  be  utilized. 

This  chapter  examines  the  mathematics  of  determining  the 
probability  density  function  of  a  performance  measure,  first 
for  performance  measures  which  are  explicit  functions  of  the 
component  variables,  then  for  implicit  performance  measures. 
Two  integration  schemes  are  then  presented  for  performing  the 
integration  on  a  digital  computer. 

A.   EXPLICIT  CASE 

Consider  a  performance  measure,  Z  as  a  function  of  the 

T  T  . 

system's  component-variables,  X  ,  where  X   is  the  transpose 

of  an  n-vector  of  elements  X, ,  X„ ,  ...,  X  .   It  is  assumed 

12         n 

that  the  X. 's  are  statistically  independent,  so  that  the 

T 
joint  probability  density  function  of  X   is  given  by  the 

product  of  the  probability  density  functions  of  the  individ- 
ual X.  's  . 

3. 


T 
The  probability  density  function  of  Z  =  g(X  )  can  be 

written  as  [7] 


pz(z)     = 


vP^T/Z(xT,z)dV  (1) 


T 

T       (x      z ) 
where   p      ,Z    —   '         is    the    joint   probability   density    function 

T 
for   X      and    Z,    V   is    the    n-dimensional    space    determined   by    the 

component-variables,    g    is    the    functional    relation   between    Z 

T 
and   the    component    variables,    X    ,    and   dV  =   dxndx„...dx    . 
r  —  1      2  n 

If   a   transformation    of   variables    is    defined   such    that 

YT      =       (X2,X3, . ..,Xn)     ;       XT      =       (XjY1)  (2) 

X±      =      f(Z,YT)  (3) 

where    f    is    the    functional    relation   between    the    component- 
variable   X,    and   the   performance    measure    Z,    equation    (1) 


becomes    [7,    9] 


pz(z) 


pvT[f (z,yT) ,yT] | J(z) | du  (4) 

U      -  -        - 


where    U   is    the   n-1   dimensional   space    determined   by   the    n-1 

T 
component-variables    Y    ,    dU   =    dy,dy«...dy    _,  ,    y.    =   x-;  +  i/    anc^ 

J(z)     is    the   Jacobian    of    the    transformation    (2)     and    (3). 
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J(z)  can  be  written  as 


J(z)   = 


9x-   9x, 


9x 


9z 

9x, 

9~i~ 


9y- 


""  9~ 


n-1 


9x 


n 


9x 


n 


9z 
9x. 


Jn-1 
0 


9x 


n 


9x. 

¥z~ 


(5) 


I  9z 

where  I  is  the  n-1  by  n-1  identity  matrix. 

Since  the  joint  probability  density  function  of  statis- 
tically independent  random  variables  is  the  product  of  the 
individual  probability  density  functions,  equation  (4) 
becomes 


Pz(z) 


r  T         T        i  T 

p       [f(z,y    ) ,y    )       9x,/9z|    pvT(y    ) dU 

U      Xl  ~  L  -     ~ 


(6) 


In   the    case   where    the    performance    measure    is    a  known    or 
explicit    function    of   the    component-variables,    the    inverse 
functional    relation    f    relating   X,    to   Z    can   usually   be   ex- 
pressed   analytically    as    can   the    partial   derivative    of    Z   with 
respect   to   X,    or    any    of   the    other    component-variables.       Cir- 
cumstances  may    require    reordering   the    component-variables 
to   obtain   these    relations    in    a    convenient    form. 
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As    an   example,    consider   the    series    R-L-C    circuit   in 
Figure    1.      The    current    I    through   the    circuit    is 

I(*"J       "       R+    Vdl-    1/o.C)  (7» 

and   is    a   maximum   for   wL   =    1/ooC.      The    angular    frequency    at 
which   the    current    falls    to    3    db   of   the   maximum   occurs   when 

R      =      coL    -    1/coC  (8) 

or, 


w3db      =      Ry/2L   -  *    r2/4Ij2    +    1/LC  O) 

Let   the    performance   measure    Z    be    the    upper    3    db    frequency 

T 
and   X      =    (R,L,C) ,    so   that 


Z   =   g(XT)  =  X1/2X2  +  /  X12/4X22  +  1/X2X3  (10) 

and  from  ( 8) , 

X1      =   f(Z,YT)   =   ZX2  -  1/ZX3  (11) 

so  that  the  partial  derivative  is 

8X1/3Z   =   X2  +  1/Z2X3  (12) 

If  each  of  the  components  has  a  Gaussian  distribution 

about  its  mean  such  that 

—   2     2 
exp [- (x  .  -x  . )  /2a  .  ] 

Px.   -     (2^l/2„3      ]   '  3=1,2,3         (13) 
i        (2tt  )    a  . 

where  x.  is  the  expected  or  average  value  of  the  variable 

X.,  and  a.  is  the  standard  deviation  of  the  variable  X . . 
D        1  3 
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swmv — \{ 


^0)  ttico) 


>'li;«> 


R 


Figure    1.       Series    R-L-C   Circuit 
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T  T 

Substituting    for    X,    in    (13)     for    j    =    1,    p       [f(z,y    ) ,y    ] 

1  x1  -       - 

becomes 

—      2  2 

exp  [- (zx0-l/zx_.-x, )    /2a,     ] 

PX        =      W23        —  (14) 

Xl  (27T)1/Za1 

T 
Since    X?    and   X      are    independent,    p  T  (y_   )     is    the    product    of 

pv    (x„)     and   pv    (x_) ,    and    (6)    becomes 
X2      2  X3      j 

—      2  —      2  —      2 

-(zx2-l/zx3-x,)  (x2~xn)  (x3~x3) 

exp[ —2 —2 —2— ] 

f°°    f00  2al  a2  3 


PZU)   = 


dx_dx. 


—  oo      — oo 


(2Tr)3/2a1a2a3[l/(x2+l/z2x3)]  2      3 

(15) 


With  a-priori  knowledge  of  the  probability  density  func- 
tions of  the  component-variables,  (6)  or  (15)  can  be  inte- 
grated on  a  digital  computer  using  one  of  many  algorithms 
available.   The  two  quadrature  formulas  presented  below  offer 
routines  which  are  fast,  and  for  the  number  of  points  at 
which  the  integrand  is  evaluated,  two  of  the  most  accurate 
routines  available  [5,  8]. 

B.   IMPLICIT  CASE 

If  the  performance  measure  is  not  known  explicitly,  but 
is  obtained  implicitly  from  a  solution  of  a  mathematical 
model  of  the  system,  relation  (3)  cannot  be  determined.   One 
may,  however,  resort  to  the  cumulative  distribution  obtained 
from  (1)  by  [7] 
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(Z     (  T 

F(z)       =      P[Z    <_  z]    =  PxT,z(x    ,z)dVdz  (16) 

— oo 


Since 


PxT,z(xT/Z)       =      pz/xT(z/xT)pxT(xT)  (17) 

T 
where   p     ,  T(z/x    )    is    the    conditional   probability   of    Z    given 

T  T  T 

that   X      has    occurred;    and    for    a    given    x    ,    p     ,  T(z/x    )     is 

either    0    or   1,    the    cumulative    distribution    can   be    approxi- 
mated by 


Z<z 


F(z)       =         I    Az  p    T(xT)dV  (181 


V      - 

where  the  integration  over  the  space  V  may  be  performed 

utilizing  the  quadrature  formulas  below  with  the  value  of 

T  T  T 

p  T  (x  )  as  calculated  for  Z  (x  )  <_  z  and  p  T(x  )  =  0  if 

A   —  —  A    — 

Z(xT)  >  z. 

The  cumulative  distribution  function  for  the  performance 
measure  can  then  be  fitted  to  a  polynomial,  and  if  desired, 
the  polynomial  can  be  differentiated  analytically  to  obtain 
the  probability  density  function. 

C.   QUADRATURE  SUMS  FOR  EVALUATING  INTEGRALS 

The  form  of  the  performance  measure  has  been  reduced  to 
an  integral  or  multiple  integral  which,  in  general,  cannot  be 
evaluated  in  closed  form.   All  algorithms  which  approximate 
the  value  of  an  integral  by  a  linear  combination  of  values  of 
the  integrand  are  exact  [5]  for  the  integrand  being  of  a 
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certain  degree  or  less;  and  in  most  cases  the  degree  is  one 
less  then  the  number  of  points  at  which  the  integrand  is 
evaluated.   Quadrature  sum  algorithms  utilize  orthoganal 
polynomials  to  arrive  at  the  form  of  the  linear  combination 
of  values  of  the  integrand;  and  as  a  result  are  exact  for  the 
integrand  being  a  polynomial  of  degree  one  less  than  twice 
the  number  of  points  at  which  the  integrand  is  evaluated. 
[5] 

Consider  the  integral 

rb 


i  = 


h(x)g(x)dx  (19) 

a 


where  a  and  b  are  real  numbers,  finite  or  infinite,  g(x)  is 
an  arbitrary  function  and  h(x)  is  a  particular  weighting 
function  to  be  described.   Then  by  selecting  a  polynomial 
Q  (x)  of  degree  n  such  that  h(x)Q  (x)  is  orthoganal  to  x  , 
for  all  m  =  0,  1,  ...,  n-1,  so  that 

x^UJQ  (x)dx   =   0  (20) 

J  a 

If  (20)  is  approximated  by  a  linear  combination  of  the 

integrand 

rb  n 

xinh(x)Qn(x)dx   =    I    A.  x,mQ(x,  )  (21) 

J  a  k=l    K 

then  the  right  hand  side  of  (21)  will  be  equal  to  zero  for 

any  set  of  values  of  A,  if  the  x,  '  s  are  chosen  such  that  they 

are  the  zeros  of  Q  (x) . 

n 
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By  forming  the  n  sums 

rb  n 

xmh(x)dx   =    I    A^™  (22) 

*  a  k=l 

for  m  =  0,1, ; . . ,n-l,  and  for  the  n  values  of  x,  being  the  n 

zeros  of  Q  (x) ,  the  n  values  of  A,  can  be  found  such  that 

rb  n 

h(x)g(x)dx   =    I    \g(x.)  (23) 

J  a  k=l 

is  exact  for  all  polynomials  of  degree  2n-l.  [5] 

It  can  be  shown  that  if  the  polynomial  Q  (x)  is  known 

j  n 

then    the    values    of   A,     can    be    found    to  be     [5,    8] . 

Ak      =      2/(l-xk2) [dQn(xk)/dx]2  (24) 

1.       Gauss    Legendre    Quadrature 

The    first    obvious    choice    of   w(x)     is    1,    and   if   the 

interval    la,    b]     can   be   normalized   to   the    interval    [-1,    1]    by 

a   change    of   variables,     (2  3)    becomes 

rl  n 

g(x)dx      =         I    Akg(x    )  (25) 

and  the  polynomial  Q  (x)  is  a  Legendre  polynomial  of  degree  n 


[5]  . 


n  ,  2  ,  v  n 


(x)     =    cr(x--i) 


Q(x)   =     ;  „  '  (26) 

n         ,  n~n  , 
dx  2  nl 

The  n  values  of  x,  can  be  found  from  the  zeros  of 

k 

Q  (x)  and  the  values  A,  from  (2  4)  or  the  values  can  be  ob- 
tained from  Table  II  for  values  of  n  from  2  to  6  and  from 
[5]  for  values  of  n  from  2  to  48.   The  magnitude  of  the  error 
bound  can  be  shown  to  be  15], 


17 


02n  +  l ,  .  >  2      I df   i 
en    <-  2     (nl)    max  — y— 
GL  <  Tt — j_ i  \  /-» — \"~r      i  2n'  (2  7) 

—  (2n+l)  (2n)  I   x   dx  v   ' 

2 .   Gauss  Hermite  Quadrature 

As  is  often  the  case  when  dealing  with  probability 

density  functions,  a  variable  is  Gaussian  in  distribution, 

and  as  a  result  an  integral  of  the  form 

f°°  n 

h(x)exp(-x2)dx   =    \    A,h(x,  )  (28) 

J  k=l  K    K 

— oo 

must  be  evaluated.   This  is  the  form  of  (19)  where  w(x)  = 
exp (-x  )  . 

The  class  of  Chebychev-Hermite  polynomials  are 

2 
orthogonal  with  respect  to  exp(-x  )  over  the  real  line,  and 

are  defined  by  [5] 

Q  (x)   =   (-l)n  exp(x2)^—  [exp(-x2)  ]      (29) 

dx 

The  values  of  x,  can  be  found  from  the  zero's  of 

k 

(29)  and  those  of  A,  from  (24) .  However  these  values  have 
been  tabulated  and  are  included  in  Table  I  for  values  of  n 
from  2  to  6  and  in  [5]  for  values  of  n  from  2  to  32. 

The  magnitude  of  the  error  bound  for  approximating 
(28)  by  a  Gauss  Hermite  Quadrature  sum  can  be  shown  to  be 
[5] 

n!  (tt)  '   max  |  d   h(x)  .  /OA* 

e   = — [ — ^ — !: — -|  (30) 

2n(2n)!    X  dxZn 
3 .   Application  to  the  Example 

The  integrand  for  the  example,  equation  (15)  is  of 
the  Gauss  Hermite  form  since  the  variables  all  have  Gaussian 
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TABLE    I 

GAUSS    HERMITE    QUADRATURE 
COEFFICIENTS    AND    ZEROS 


Number 

of  Evaluation 

Zeros 

;    x  . 

Coefficients    A. 

Points , 

n 

D 

3 

2 

+0.70  7 

1068 

0.886    2269 

3 

+1.224 
0.0 

745 

0.295    4090 
1.181    635 

4 

+1.650 

680 

0.081    31284 

+0.524 

6476 

0.804    9141 

5  +2.020  183  0.019  95324 
+0.958  5725  0.393  6193 

6  +2.350  605  0.004  530010 
+1.335  849  0.157  0673 
+0.436  0774  0.724  6296 
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TABLE    II 

GAUSS    LEGENDRE    QUADRATURE 
COEFFICIENTS    AND    ZEROS 


Number 

of  Evaluation. 

Zeros 

;    x  . 

Coefficients    Au 

Points , 

n 

3 

J 

2 

+0.577 

3503 

1.000 

0000 

3 

+0.774 

5967 

0.555 

5556 

~0.0 

0.888 

8889 

4 

+0.861 

1363 

0.347 

8548 

+0.339 

9810 

0.652 

1452 

5 

+0.906 

1798 

0.236 

9269 

+0.538 

4693 

0.478 

6287 

0.0 

0.568 

8889 

6 

+0.932 

4695 

0.171 

3245 

+0.661 

2094 

0.360 

76  16 

+0.238 

6192 

0.467 

9139 

20 


distributions.   By  changing  variables,  let 


x. 


x2  X2 

/To  n 


(31) 


x. 


x3  X3 

/To  ^ 


(32) 


dx2 


dx, 


/To. 


(33) 


dx' 


dx. 


/2a. 


(34) 


Since  X2  and  X..  are  independent,  the  double  integral 
becomes  the  product  of  the  two  integrals  and  (15)  takes  the 
form 


oo  ,0O 


I  2         2 

g(z,x^,x^)exp(-xj    -x^    ) dx^dx* 


— oo         _oo 


n        n 


I         I    A.Akg(z,x       ,x      ) 
j=l  k=l    J  j         k 


(35) 


where 


g(z,x^,x3) 

exp[-{z(2"L/     ?2x^+x2) -l/z(2x/  'a0x^+x3) -x-j^}^  6] 


1/2,,       ._,_-  ,     w„Ml/2^       ,.-■,  ,2,^2 


-1/2,    ,3/2       nl/2  ,    -  ^,  /   2  ,~l/2  ,    -   ,  ,  -1 

2    /     (tt)    '    a,  [2    '    a2x'+x  +l/z    (2    '    a-xi+x-J  ] 


(36) 


The  value  of  n  can  then  be  chosen  as  desired  and  the  sum 

in  (35)  formed  using  the  values  of  A,  and  x„   and  x-.   from 

K       ^k       \ 

Table  II  or  [5].   The  expression  (36)  becomes  quite  formidable 
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to   differentiate    after   two   or  more    differentiations,    so   that 
one   must    assume    that   the   error   bound    (30)    is    small.      The 
approximation    (35)    has    been   performed    for   n   =   6    and   the    mean 
values    R  =    1000    ohms,    L  =    1.0    millihenry,    C  =    0.001   micro- 
farad.      The    standard   deviation    of   each    variable   was    assumed 
to   be    3.33%    of   its   mean    value.      The    results    are    contained   in 
Chapter   IV. 

Not    all   variables   have    Gaussian    distributions,    and  most 
performance   measures,    if    functions    of   many   variables,    are 
implicit    functions.       If   the   example    performance    measure    is 
considered   to   have    been    obtained   implicitly,    with    the    same 
density    functions    for   the    components,    the    performance   meas- 
ure's   cumulative    distribution    function    can   be    determined   by 
approximating    (19)    by    a   Gauss    Legendre   Quadrature    sum.      Since 
exp(-4.5)    is    very   nearly    zero,    the    limits    of    -°°    and   °°    can   be 
changed   to    -3a    to    3a    for  each    variable    and   by   the    following 
change    of   variables,    normalized   to  the    interval    [-1,    1].      Let 


x 


x.  -x. 

l      r 

3a. 

l 


i   =    1,    2,    3 


(37) 


dx.' 

r 


so  that  (18)  becomes 


F(z)   = 


dx./3a. 


fl   rl 
£  Azg(x' ,x' ,x' )  dx'dx' dx' 


-1   -1 


(38) 


(39) 


where 
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2  7exp(-x,2/2-x'2/2-x'2/2) 

g(x'  x'  x')   = =-7^ (40) 

J  (2tt)j/z 

The  approximation  of  (39)  by  the  summation 

n    n    n 
F(z)   =1111  AzA  A  -A.  g(x'  ,xl  ,X«  )  (41) 

z<z  i=l  j=l  k=l   L  D  K    1  zj      \ 

just  requires  the  n  values  of  A,  and  the  n  values  of  x!   and 

K  1k 

evaluating  g  at  these  values  of  x!  .   It  must  be  remembered 

Xk 

T 
that  g  takes  the  form  of  (40)  only  for  Z  (x  )  <,   z  for  some  z; 

otherwise  g(x')  =0. 

Equation  (41)  is  evaluated  for  n  =  6  and  the  same  values 

of  the  components  as  above  for  the  Gauss  Hermite  example,  and 

the  results  are  contained  in  Chapter  IV  for  comparison.   As  in 

the  Gauss  Hermite  example,  the  error  is  assumed  to  be  small 

due  to  the  factorials  in  the  denominator  of  (27) . 
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III.   MONTE  CARLO  METHODS 

This  chapter  examines  the  most  crude  of  the  Monte  Carlo 
techniques  for  obtaining  statistics  on  a  performance  measure, 
to  methods  which  attempt  to  approximate  single  integrals  in 
such  a  manner  that  the  approximation  is  reasonably  accurate. 
The  approximations  are  then  extended  to  apply  to  multiple 
integrals  in  such  a  manner  that  fewer  evaluations  of  the 
integrand  are  made  than  in  conventional  methods. 

Monte  Carlo  methods  consist  of  solving  problems  of  a 
computational  nature  by  constructing  a  random  process  for  the 
problem,  such  that  the  parameters  of  the  random  process  are 
the  desired  quantities  of  the  problem. 

Random  processes  usually  imply  random  variables  which 
must  be  generated  in  a  manner  as  to  represent  typical  values 
from  the  variable's  probability  density  function.   Random 
number  generation  is  then  an  important  aspect  in  any  Monte 
Carlo  method  of  approximation.   Although  random  number  gener- 
ation is  a  field  of  its  own,  Appendix  A  treats  the  subject 
suitably  for  purposes  of  this  thesis. 

A.   CRUDE  MONTE  CARLO 

The  most  crude  form  of  Monte  Carlo  is  that  of  taking  sam- 
ples of  the  performance  measure  to  obtain  an  expected  value 
and  perhaps  a  variance  of  the  performance  measure.   If  the 
component-variables  are  generated  in  such  a  way  so  that  the 
values  are  representative  of  their  distribution  functions,  and 
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for  each  set  of  values  X   of  component  variables  generated, 
the  performance  measure  Z.  =  ZCX1)  is  evaluated,  then  the 
expected  or  average  value  of  the  performance  measure  is  [7] 

(Z.)  =  1  j   I  (42) 


n 
Z   =   E 


where  n  is  the  number  of  samples. 

While  for  sufficiently  large  n,  Z"  is  an  accurate  estimate 
of  the  average  value  of  the  performance  measure,  no  knowledge 
of  the  form  of  the  probability  density  function  of  Z  is  de- 
rived; and  consequently  no  knowledge  of  the  likelihood  of 
other  values  of  Z  is  obtained. 

A  more  useful  scheme  would  be  to  order  the  sample  values 

of  Z .  as  follows  to  obtain  a  cumulative  distribution  for  the 

l 

performance  measure. 

Consider  a  random  variable  S  such  that  if  a  value  of  the 
performance  measure  Z  is  sampled,  and  the  sample  value  Z.  is 
less  than  an  arbitrary  value,  say  a,  then  S.  =  1;  and  if  Z. 

is  greater  than  the  value  a,  S.  =  0.   Then  S  represents  a 

T 
random  process.   If  the  values  of  the  component-variables  X 

are  sampled  from  their  density  function  as  above,  then  the 

cumulative  distribution  function  for  Z  can  be  approximated  by 

l   n 
F(a)   =   P[Z<a]  =  -   T  S.  (43) 

—     n  .  L,    l 
i=l 

where   n   is    the   number   of   times    the   performance    measure    is 
sampled. 

Since   each   of   the    component-variables    are    sampled   inde- 
pendently   for   each    sample    value    Z.,    the    number    of   samples    is 
not   directly   dependent   on    the   number   of   variables. 
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This  crude  method  yields  a  cumulative  distribution  func- 
tion which  may  be  fitted  to  a  polynomial  and  differentiated 
analytically  to  obtain  the  probability  density  function  for 
the  performance  measure  Z. 

As  an  example  of  the  two  methods,  consider  the  example 
of  Chapter  II.   Since  the  three  variables  X,,  X? ,  and  X-.  are 
assumed  to  be  Gaussian,  three  variables  w,  ,  w„ ,  and  w->  are 
generated  by  the  method  in  Appendix  A  to  approximate  Gaussian 
distributions,  each  with  0  mean  and  variance  of  1.   The  val- 
ues of  the  samples  x, ,  x0 ,  x_  are  then 

x.   =   w.a.  +  x.  ,   i  =  1,2,3  (44) 

l       rii  ' 

where    a.    and   x.    are    the    standard   deviation    and  expected   val- 
11  c 

ue ,  respectively  of  the  variable  w..  These  values  are  then 
used  to  evaluate  the  performance  measure.  If  the  expected  or 
average  value  of  the  performance  measure  is  desired,  the  cal- 
culated values  of  the  performance  are  averaged  as  in  (42) . 
If  a  cumulative  distribution  is  desired,  then  for  k  values  of 
a.,  j=l,2,...,k  the  random  variable  S  defined  above  is  summed 
to    form   F(a.) .       These    values    for   the   examole    are    listed   in 

y 

Chapter  IV,  and  a  computer  program  listing  for  the  computa- 
tion is  provided  in  Appendix  B. 

B.   CONFIDENCE  INTERVALS 

The  crude  Monte  Carlo  method  above,  as  with  all  Monte 
Carlo  methods,  yields  results  which  are  as  accurate  as  de- 
sired within  a  probability  or  confidence  determined  by  the 
desired  accuracy,  and  the  method  of  approximation  used.   What 
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follows  is  a  statistical  estimate  of  the  error  introduced  by 

use  of  Monte  Carlo  approximations. 

If  A  is  an  event  (such  as  the  event  that  the  values  of  a 

performance  measure  is  less  than  some  specified  value) ,  and  R 

is  a  random  variable,  such  that  R.  =  1  if  A  occurs  on  the  ith 

sample,  and  R.  =  0  otherwise,  and  if  T  is  a  random  variable 

such  that 

n 
T   =    I      R.  (45) 

i=l   x 

for  n  samples,  then  T  is  the  number  of  times  that  the  event  A 

occurs  in  n  samples.   The  expected  value  of  T  is 

n 
E(T)   =   E(  I    R.)  =  nE(R.  )  =  np  (46) 

i=l  *■ 

where  p  is  the  probability  of  the  event  A  occurring.   The 

variance  of  T  is 

n        9  9 
V(T)   =   V(  J  R.)  =  no  (47) 

i=l  x 

2 

where  a   is  the  variance  of  the  event  A. 

From  Chebychev's  inequality  [7], 

P      =      P[|T-np|<    d]     >    1    -    a2/nd2  (48) 

where  d  is  an  arbitrarily  small  number.  Equation  (48)  states 
that  for  a  fixed  confidence  of  |T-np|  being  less  than  a  small 
value    d,    the   number   of    samples    to    achieve    the    accuracy,    d, 

is    bounded  by 

2 

n    >    -^ (49) 

~  dZ(l-P) 

where  P  is  the  desired  confidence  or  probability. 
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Another  way  of  looking  at  (4  8)  is  that  for  a  fixed  con- 

-1/2 
fidence  P,  the  error  of  the  estimation  varies  as  an     .   If 

the  variance  of  the  random  variable  can  be  reduced  to  0,  the 

confidence  of  the  estimation  could  be  1  for  an  error  of  0. 

For  crude  Monte  Carlo,  the  random  variable  S,  which  takes 

on  values  of  1  whenever  the  sampled  performance  measure  is 

less  then  a  specified  value,  and  values  of  0  otherwise,  has  a 

binomial  distribution.   If  p  is  the  expected  value  of  S,  then 

2 

the  variance  of  S  is  a   =  p(l-p)  [7] .   Thus  for  a  confidence 

of  P  =  0.9  of  being  within  d=0.01  of  np  for  crude  Monte  Carlo, 

n  >~    (oioo'l)  (0.1)   ■   P<1-P>1°5  <50) 

Since  the  maximum  value  of  p(l-p)  is  0.25,  inequality  (50) 
states  that  n  must  be  greater  than  twenty-five  thousand  in 
order  to  insure  the  accuracy  of  0.01  with  a  confidence  of 
0.9. 

C.   VARIANCE  REDUCTION  BY  ESTIMATING  INTEGRALS 

As  was  seen  in  Chapter  I,  the  difficulty  of  determining  a 
performance  measure's  probability  density  function  lies  in 
evaluating  the  integral  (1).   Two  routines  were  presented 
which  approximated  the  integral  with  some  accuracy.   The 
crude  Monte  Carlo  method  illustrated  above  can  be  thought  of 
as  an  approximation  to  an  integral.   Other  methods  exist, 
however  which  reduce  the  variance  by  hundreds  of  thousands 
over  crude  Monte  Carlo  [3] .   These  methods  are  presented  below 
for  single  integrals  and  then  a  general  extension  is  made  for 
multiple  integrals. 
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If  I  is  the  estimator  for  the  integral, 


f  (x)  dx,  where 
0 


x  is  a  variable  normalized  to  the  interval  [0,  1] ,  then  by 
letting 

I   =  ±{^\  (51) 

py(y) 

where  pY(y)  is  the  probability  density  function  of  the  ran- 
dom variable  Y  and  f(y)  is  the  value  of  the  function  f  at 
some  point  y  sampled  from  pY(y)/  the  expected  value  of  I  is 

E(I)   =  E[^{VA  ]   =   f   f(y)dy  (52) 

It  is  then  necessary  to  select  the  probability  density  func- 
tion pY(y)  in  such  a  manner  as  to  minimize  the  variance. 

Define  a  function  f*(x)  such  that  the  function  which  is 
f(x)-f*(x)  is  a  straight  line  passing  through  the  end- 
points  of  f (x)  on  the  interval  [0,  1] ,  as  shown  in  Figure  2. 
That  is, 

f*(x)   =   f(x)  -  (l-x)f(O)  -  xf(l)        (53) 

The  difference  function  f(x)-f*(x)  can  be  readily  inte- 
grated and  can  be  used  as  a  first  approximation  to 

fl 

f (x) dx  as 

0 

1 

[f (x)-f*(x) ]dx   =   l/2f(0)  +  1/2 f(l)    (54) 
0 


Now    if    I*    is    defined    as 

I*      =      *±&)  (55) 

pY(y) 

it  becomes  necessary  to  sample  f*(y)  to  get  the  estimate  from 
(53)  ,  (54)  and  (55)  for  I 
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Figure    2.       Difference    Curve   Estimation. 
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T       f (0)  +  f (1)     f (y)  -  (l-y)f(Q)  -  y(l) 

2        +         pY(y)  ,       <56> 

To  determine  py(y),  let  f(x)  be  quadratic  in  x,  that  is 

f(x)   =   Ax2  +  Bx  +  C  (57) 

By  substituting  for  f  in  (56)  the  f  of  (57) ,  and  replacing  I 
with  the  integral  of  f  in  (5  7) ,  that  is 

I   =   A/3  +  B/2  +  C  (58) 

one  can  obtain  for  pY(y) , 

PY(y)   =   6y(l-y)  (59) 

It  can  be  shown  by  substitution  that  the  approximation  I  in 
(56),  with  pY(y)  of  (59),  is  a  zero  variance  estimator  for 
the  integral  of  quadratic  integrands. 

The  function  f * (x)  could  have  been  defined  as 

f*(y)   =   f(l-y)  -  yf(0)  -  (l-y)f(l)  (60) 

with  a  resulting  estimator  of 

t      f(0)  +  f(D  .  f(l-y)  -  yf(0)  -  (l-y)f(l) 
2  2  6y(l-y) 

(61) 

which  is  also  exact  for  f (x)  quadratic  in  x.   Since  y  and 
1-y  have  the  same  distribution,  but  lie  on  the  opposite  side 
of  1/2,  the  average  of  (56)  and  (61) ,  as  can  be  shown  by 
substitution,  results  in  a  zero  variance  estimator  for  f(x) 
cubic  in  x[3] .   Thus 

t       f(0)  +  f(l)  ,  f(y)  +  f(l-y)  -  f(0)  -  f(l)     (ao. 
Z3  ~      2 +  12y(l-y) (62) 

I,  is  exact  for  f(x)  cubic  in  x. 
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In   the    same   manner   as    above,    by   defining    f*    as 

f*(x)      =      f  (x)-(l-x)  (l-2x)f  (0)+x(l-2x)  f(l)-4x(l-x)  f  (j) 

(63) 

a  difference  curve  which  is  quadratic  and  passes  through  the 
two  end  points  and  the  center  of  f(x)  becomes  the  basis  for 
the  initial  estimate  for  the  integral.   Applying  the  methods 
above,  a  zero  variance  estimator  for  quartic  integrands  be- 
comes { 3]  . 

t       f (0)+f (l)+4f (1/2) 
X4  6 

+  f  (y)-d-y)  (l-2y)  f  (0)  +y  (l-2y)  f  (1)  -4y  (1-y)  f  (1/2) 

30y(l-y) (l-2y)2 

(64) 

and   a   zero   variance   quintic   estimator   becomes    2 

t         =      f  (0?+f  (D+4f  (1/2) 
5  6 

+   f  (y)+f  (i-y)-d-2y)2  [f  (0)+f  d)1  -8y  ( 1-y)  f  (1/2) 

30y(l-y) (l-2y)2 

(65) 

The  probability  density  function  for  both  the  quadratic 
and  cubic  zero  variance  estimators  is  (59)  while  that  for  the 
quartic  and  quintic  zero  variance  estimators  is 

py(y)   =   30y(l-y) (l-2y)2  (66) 

D.   EXTENSIONS  TO  MULTIPLE  INTEGRALS 

The  zero  variance  estimators  above  were  derived  in  [2] 
only  for  single  integrals.   While  the  method  can  be  extended 
to  multiple  integrals  by  nesting,  nothing  is  gained  since  the 
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f(x,y)dxdy,  consider  the 
0 


number  of  function  evaluations  increases  exponentially  with 
the  number  of  variables  of  integration.   In  this  section, 
double  integral  extensions  to  the  quadratic,  cubic,  quartic, 
and  quintic  zero-variance  estimators,  and  triple  integral 
extensions  to  the  quadratic  and  cubic  estimators  are  proposed, 
and  the  method  of  derivation  is  presented.   Although  not 
thoroughly  tested  on  a  wide  range  of  integrands,  it  is  be- 
lieved that  this  method  of  extension  can,  for  multiple  inte- 
grals with  many  variables  of  integration,  provide  an  accurate 
means  of  approximating  the  integral  with  fewer  evaluations  of 
the  integrand  than  for  the  quadrature  sum  routines. 

f1  I 
For  the  double  integral 

J0 
difference  function  which  is  a  surface  passing  through  the 
four  corners  f(0,  0),  f(0,  1),  f(l,  0),  and  f(l,  1).   This 
function  is 

f (x,y)-f*(x,y)   =   (1-x)  (l-y)f(0,0)  +  (l-x)yf (0,1) 

+  x(l-y)f(l,0)  +  xyf(l,l)        (67) 

Integrating  this  difference  function  results  in 

f1  f1  r*t         ^  **/    \ia    a                f  (0,0)+f  (1.0)+f  (0,l)+f  (1,1) 
[f  (x,y) -f*(x,y)  ]dxdy   =   '- 2 

J0    J0 

(68) 

Now   if    I*    is    defined   similar   to    (55)     as 

T*      =      f*(x*y> (69) 

Px(x)py(y) 

and  pv(x)  and  pv(y)  have  the  form  of  (59),  solving  for 


33 


f*(x,y)  in  (67)  and  using  (69)  as  the  estimate  for 
f 


0 
1 


-  f * (x,y) dxdy ,  the  estimate  to        f(x,y)dxdy  becomes 
0  J0 


0 


f  I      *,         ,,    ,                f  (0,0)+f(0,l)+f (l,0)+f (1,1) 
f(x,y)dxdy   =   — - — - — — - — — — i — - — - — - — - 

J0  J0 

,  f (x,y)-(l-x) (1-y) f (0,0)-(l-x)yf (0 f 1) -x(l-y) f (1,0) -xyf (1,1) 

36xy(l-x) (1-y) 

(70) 

which  is  the  double  integral  extension  to  (56) .   Replacing  x 
with  1-x  and  y  with  1-y  in  (70)  and  averaging  the  resulting 
expression  with  (70)  gives  the  double  integral  extension  to 
the  cubic  zero  variance  estimator  (62) 

1  I1    f(x,y)dxdy   =   f(0,0)+f(0,lKf(l,0)+f(l,l) 


0 


0 


f  (x,y)+f  (l-x,l-y)-(l-x-y+2xy)  [f(0,0)+f(l,D] 

72xy(l-x) (1-y) 

-(x+y-2xy) [f (0 , 1) +f ( 1, 0) ]  , 

72xy(l-x) (1-y)  K  'L) 

For  higher  order  extensions,  a  difference  surface  is 

defined  such  that  it  passes  through  the  2   end  points  of 

T  T 

f(x  ),  where  x   is  a  k  vector.   This  difference  curve,  which 

is  a  polynomial  in  each  of  its  variables,  can  be  integrated 

analytically  to  give  a  first  approximation  to  the  desired 

integral.   The  estimate  is  then  refined  by  sampling  the  f* 

function  weighted  by  the  distribution  of  the  variables,  as  in 

(69).  Thus  the  triple  extension  to  the  quadratic  estimator 

(56)  becomes 
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f1  f1  f1  f(x,y,z)dxdydz  =  f(0*0,0)+f(0,0,l)+f(0,l,Q)+f(0,l,l) 

Jo  Jo  Jo  8 

+  f(l,Q,0)+f(l,0,l)+f(l,l,0)+f(l,l,l)         f (x,y,z) 

8  216xyz (1-x) (1-y) (1-z) 

.  -d-x)  (1-y)  (1-z)  £(0,0,0)- (1-x)  (1-v)  zf  (0,0,1) 
216xyz (1-x) (1-y) (1-z) 

-(l-x)y(l-z)f (0,1,0) 
216xyz (1-x) (1-y) (1-z) 

-(l-x)yzf (0,1,1) -x( 1-y) (1-z) f ( 1, 0 , 0 ) -x (1-y) zf ( 1, 0 , 1) 

216xyz (1-x) (1-y) (1-z) 

+  -xyd-z)  f  ( 1,1,0  )-xyzf  ( 1,1,1)  ,79, 

216xyz(l-x) (1-y) (1-z)  t  { '*' 

and  the  triple  extension  to  the  cubic  estimator  (62)  becomes 
'X  f1  *f„   ..  ^ *„*„*„    _  f(0,0,0)+f(0,0,l)+f(0,l,0)+f(0,l,l) 


r1 


o 


f (x,y , z) dxdydz  = 


0  ^0  8 


,  f (1, 0,0) +f( 1,0,1) +f( 1,1,0 )+f( 1,1,1)  f(x,y,z)+f (1-x, 1-y, 1-z) 

8  432xyz (1-x) (1-y) (1-z) 

-(1-x-v-z+xz+yz)  [f(0,0,0)+f(l,l,l)] 
432xyz (1-x) (1-y) (1-z) 

-(z-xz-yz+xv) [f(0,0,l)+f(l,l,0)] 
432xyz (1-x) (1-y) (1-z) 

-(y-xy-zy+xz) [ f (0 , 1 , 0 ) +f ( 1 ,0 , 1) ] 
432xyz (1-x) (1-y) (1-z) 

-(x-xz-xy+yz)  [f(l,0,0)+f(0,l,D]  n^ 

432xyz(l-x) (1-y) (1-z)  K       ' 

For  the  quartic  and  quintic  double  integral  extensions, 
the  difference  surface  must  additionally  pass  through  the 
four  points  which  are  peripheral  midpoints  and  the  internal 
midpoint  of  the  integrand.   The  double  integral  extension 
for  the  quartic  zero  variance  estimator  becomes 
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f(0,0)+f(0,l)+f(l,0)+f(l,l) 

f  (x,y)  dxdy   = 


0  36 

f(oi)+f(ii)+fd,  o)+f  (i,i)+4f(y,y) 

+  4  [ ± ± £ 1 z    A    ] 

36 

f(x,y)-(l-x) (l-2x) [(1-y) (l-2y) f(0,0) 

+  2 2 

900xy(l-x) (1-y) (l-2x) z ( l-2y) ^ 

+  4y(l-y)  f(0,j)-y(l-2y)  f(0,l) ] 
900xy(l-x) (1-y) ( l-2x) 2 ( l-2y) 2 

-4x(l-x)  [(1-y)  (l-2y) f (|,0) +4y (1-y) f ( \,  \)  -y (l-2y) f|,l) ] 

900xy(l-x) (1-y) ( l-2x) 2 ( l-2y) 2 

x(l-2x)  [(1-y)  (l-2y)f(l,0)+4y(l-y)f(l,y)-y(l-2y)f(l,l)  ] 

+ — (74 

900xy(l-x)  (1-y)  (l-2x) z ( l-2y ) z 

and  the  double  integral  extension  for  the  quintic  zero 

variance  estimator  becomes 

\\    \\    f(x,y)dxdy   -   f(0,0)+f(0,l>+f(l,0>+f(l-l> 

f  (0,y)+f  (l,y)+f  (y,0)+f  (y,l)+4f  (y,y) 

+  4[ £ £ £ 1 £_±_] 

36 

f  (x,y)+f  (l-x,l-y)  -(l-2x)  (l-2y)  (l-x-y+2xy)  [f(0,0)+f(l,D] 


+ 


1800xy(l-x) (1-y) (l-2x) 2 (l-2y) 2 

(l-2x) (l-2y) (x+y-2xy) [ f (1, 0) +f (0 , 1) ] 

1800xy(l-x) (1-y) ( l-2x) 2 ( l-2y) 2 
-4y(l-2x)  (1-y)  [f (0 ,  j)  -f (1,|)  ] 

1800xy(l-x) (1-y) (l-2x) 2 (l-2y) 2 


-4x(l-x)  (l-2y)  [f  (pO)-f  (jrl)  ]-32xy(l-x)  (1-y)  f(j,j) 
1800xy(l-x)  (1-y)  (l-2x) 2  (l-2y) 2 


(75) 
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For  higher  order  extensions  to  the  quartic  and  quintic 

estimators,  the  difference  function  must  be  defined  so  as  to 

k  k-1 

pass  through  the  2   corners,  the  k2     peripheral  midpoints, 

T  T 

and  the  central  midpoint  of  f  (x  )  ,  where  x   is  a  k  vector. 

These  extensions  become  difficult  to  derive  due  to  the  large 

number  of  points  involved. 

The  triple  integral  extension  to  the  cubic  estimator,  (73! 

has  been  applied  to  the  example  of  Chapter  II.   The  results 

are  contained  in  Chapter  IV  for  comparison,  and  a  printout 

of  the  program  is  contained  in  Appendix  B. 
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IV.   CONCLUSIONS 


Chapter  II  presents  the  mathematics  for  obtaining  a 
performance  measure's  probability  density  function  for  ex- 
plicit and  implicit  performance  measures.   Two  integration 
schemes  using  Gauss  quadrature  sums  were  then  discussed. 
Chapter  III  illustrated  a  crude  Monte  Carlo  method  for  ob- 
taining a  performance  measure's  distribution  and  then  dis- 
cussed the  confidence  of  such  schemes.   Monte  Carlo  schemes 
for  evaluating  single  integrals  were  discussed.   Double  and 
triple  integral  extensions  to  these  methods  were  then  derived. 

From  observation  of  the  problem,  it  seems  intuitive  that 
most  system's  performance  measures  will  be  of  the  implicit 
type,  and  in  general  the  performance  measure  will  be  a  func- 
tion of  several  variables.   In  the  process  of  evaluating  the 
performance  measure's  distribution,  the  integral  (1)  must  be 
evaluated.   Evaluation  by  the  quadrature  sum  methods,  (25) 
and  (29)  require  some  n   evaluations  of  the  integrand,  where 
n  is  the  number  of  points  in  the  quadrature  sum,  and  k  is  the 
number  of  variables.   This  number  can  become  quite  large,  and 
since  the  computation  time  for  such  routines  is  roughly  pro- 
portional to  the  number  of  times  the  integrand  is  evaluated, 
this  process  could  take  a  good  deal  of  time.   However,  this 
method  is  quite  accurate,  and  it  is  easily  programmed  for  a 
large  number  of  variables. 
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The  reduced  variance  estimators  reduce  the  number  of 
evaluations  of  the  integrand  required.   For  a  single  esti- 

mate,  the  quintic  extension  requires  evaluations  at  the  2 

k-1 
corners,  the  k2     peripheral  midpoints  and  the  3  internal 

k-1 
points,  or  a  total  of  3  +(k+2)2     points,  where  again  k  is 

the  number  of  variables.   For  additional  estimates,  only  two 

additional  points  per  estimate  are  required.   Thus  M  estimates 

k-1 
of  the  integrand  would  require  2M  +  1  +  (k+2)2     evaluations 

of  the  integrand.   This  can,  in  most  cases  be  much  less  than 
that  required  for  the  quadrature  sums. 

The  drawback  to  the  reduced  variance  scheme  is  that  the 
geometric  form  of  the  integrand  must  be  examined  to  insure 
that  the  major  portion  of  the  surface  will  not  be  neglected  in 
the  sampling  process.  As  illustrated  in  Figure  3,  the  dif- 
ference curve  does  not  provide  a  good  representation  of  the 
function,  while  splitting  the  curve  into  two  sections  as  in 
Figure  4  enables  one  to  obtain  a  good  representation  of  the 
function  by  the  use  of  two  difference  curves. 

The  generation  of  the  random  numbers  with  the  probability 

2 
density  function  30x(l-x) (l-2x)   can  be  bothersome  to  program 

as  can  the  quintic  estimator  extension.   The  cubic  estimator 
is  much  more  easy  to  extend,  and  the  density  function  6x(l-x) 
is  more  readily  programmable.   However  more  points  are  re- 
quired by  the  cubic  to  ensure  as  good  accuracy  as  the  quintic 
extension. 

Properly  used,  the  Monte  Carlo  method  of  integral  estima- 
tion can  provide  a  more  rapid  method  of  determining  a 
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1  X 


Figure  3.   Improper  Difference  Curve  Choice 
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X:o 


r,      «'*' 


Figure    4.      Proper    Difference    Curve    Choice 
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performance   measure's    density    function  with    reasonable    ac- 
curacy  than   other   methods    which    require   more    standard  inte- 
gration   routines. 

A.  RECOMMENDATIONS  FOR  FURTHER  STUDY 

Monte  Carlo  schemes  seem  to  be  attractive  for  reducing 
the  amount  of  time  required  in  the  evaluation  of  the  multi- 
ple integrals  required  in  evaluation  of  a  performance 
measure's  distribution.   Higher  order  zero  variance  estima- 
tors could  be  generated  by  fitting  a  difference  curve  to  the 
points  which  are  the  zero's  of  a  Legendre  polynomial  and 
evaluating  the  required  probability  density  function  to 
achieve  a  zero  variance  estimator  for  integrands  of  order 
2n-l,  where  n  is  the  order  of  the  Legendre  polynomial.   Other 
schemes  of  generating  random  numbers  of  particular  distribu- 
tions could  as  well  be  investigated. 

B.  RESULTS    OF    APPLICATIONS    TO    EXAMPLE 

The  cumulative  distribution  function  of  the  upper  3  db 
frequency  of  the  series  R-L-C  circuit  of  Figure  1  has  been 
obtained  by  four  methods. 

First,  using  the  analytical  expressions  (11)  and  (12)  for 
the  integral  (6) ,  the  probability  density  function,  equation 
(15)  was  obtained  by  evaluating  the  integral  with  a  6-point 
Gauss  Hermite  quadrature  sum  routine.   The  cumulative  dis- 
tribution function  was  then  obtained  from  the  probability 
density  function  by  rectangular  integration.   The  resulting 
distribution  has  been  plotted  in  Figure  5,  and  will  serve  as 
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the  reference  since  integration  of  explicit  functions  is 
generally  more  accurate  than  integration  of  implicit 
functions . 

Using  the  same  performance  measure  as  for  the  explicit 
case,  but  stipulating  that  the  value  of  the  performance  meas- 
ure was  to  be  obtained  implicitly,  integral  (18)  was  evalu- 
ated using  6 -point  Gauss  Legendre  quadrature  sums.   As  can  be 
seen  in  the  plot  of  the  cumulative  distribution  function  for 
this  case,  Figure  6,  the  distribution  is  very  nearly  the  same 
as  that  obtained  in  the  explicit  case. 

Then  for  the  crude  Monte  Carlo  case,  1000  samples  of  the 
performance  measure  (10)  were  taken  with  the  values  of  the 
three  variables  being  obtained  from  a  Gaussian  random  number 
generator.   As  can  be  seen  in  Figure  7,  the  resulting  cumula- 
tive distribution  function  is  very  nearly  the  same  as  that 
obtained  in  the  explicit  case. 

The  extended  cubic  zero  variance  estimator  equation  (73) 
was  then  applied  to  the  integral  (18) .   The  required  distri- 
butions of  the  components  were  generated  as  explained  in 
Appendix  A.   Twenty  approximations  for  each  value  of  the  per- 
formance measure  averaged  to  obtain  the  cumulative  distribu- 
tion plotted  in  Figure  8.   Two  discrepancies  are  observed 
here.   First,  the  cumulative  distribution  exceeds  the  maximum 
value  of  1,  and  secondly,  there  is  an  apparent  discontinuity 
in  the  curve.   Since  the  entire  curve  has  a  higher  value  than 
that  of  Figure  5,  it  is  assumed  that  the  difference  surface 
estimation  was  too  large. 
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Although  the  example  did  not  prove  the  worth  of  the 
Monte  Carlo  method,  a  count  of  the  required  number  of  evalu- 
ations for  high  multiple  integrands  for  a  6 -point  Gauss 
Legendre  routine  as  compared  to  that  required  in  the  extended 
cubic  estimator  averaging  10  0  approximations  for  each  value 
of  the  performance  measure  will  illustrate  its  merits.   For 
a  multiple  integral  with  10  variables  of  integration,  the 
6-point  Gauss  Legendre  method  would  require  6    evaluations 
of  the  integrand.   For  the  same  integrand,  the  cubic  exten- 
sion would  require  200  +  2    evaluations.   The  savings  in 
computation  time  should  be  apparent. 
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Appendix  A 
RANDOM  NUMBER  GENERATION 

A.   UNIFORM  RANDOM  NUMBERS 

A  truly  uniform  random  number  cannot  be  generated  on  a 
digital  computer.   However,  methods  such  as  the  power  residue 
[4,  10]  method  can  generate  a  sequence  of  numbers  which  have 
the  properties  of  uniformly  distributed  random  numbers. 

The  method  of  power  residues  is  based  on  modular  arith- 
metic.  If  X.  is  an  arbitrary  integer  in  the  sequence  of 
pseudo-random  numbers  to  be  generated,  the  next  number  gener- 
ated is  determined  by 

X.,,   =   X.p  (modulo  M)  A-l) 

i+I      r 

2 
where  p  is  any  prime  integer  such  that  p   is  greater  than  M, 

the  modular  base  of  the  computer's  arithmetic  unit.   The 

number  X./M  is  then  an  approximation  to  a  random  variable 

with  a  uniform  distribution  on  the  interval  (0,1) .   Utilizing 

a  digital  computer  with  a  32  bit  register  (31  number  bits 

plus  one  sign  bit) ,  a  value  of  p  of  65,539  for  an  M  of 

31  29 

2,147,483,647  =  2   -1,  a  sequence  of  2    numbers  can  be  gener- 
ated before  any  number  of  the  sequence  is  repeated  [4]. 

Since  the  digital  computer  automatically  performs  modular 
arithmetic  on  integers,  with  the  modulus  being  determined  by 
the  size  of  the  registers  in  the  arithmetic  unit,  uniform 
pseudo  random  numbers  can  be  generated  quite  rapidly. 
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B.  GAUSSIAN    RANDOM   NUMBERS 

Gaussian    random  numbers    are   most    commonly    generated  uti- 
lizing  uniformly    distributed   pseudo-random  numbers.       If 
R. ,    i=l,2,...,n   are    independent    samples    from   a    density    func- 
tion   uniformly   distributed   over   the    interval    (0,1) ,    then    a 
normally   distributed   random  variable   with   0    mean    and   vari- 
ance   of    12/n    can    be    approximated   by    [3,     4,    10] 

1?      n 
Z      =      —      I     (R.)     -    6  (A-2) 

n    .  *•,      i 
i=l 

If  n  =  12,  Z  will  have  0  mean  and  variance  of  1. 

C.  6x(l-x)  DISTRIBUTION 

If  R.  is  an  independent  sample  from  a  random  variable 
uniformly  distributed  on  the  interval  (0,1),  for  i  =  1,2,3, 
and  if  the  three  samples  are  rearranged  in  order  of  magni- 
tude such  that 

R-L  <  R2  <  R3  (A-3) 

then  R  has  probability  density  function  6x(l-x)  [3]. 

D.  30x(l-x) (l-2x)2    DISTRIBUTION 

If    R    ,    R    ,    R_,    R    ,    and    R      are    independent   samples    from    a 
random  variable    uniformly    distributed   over   the   interval    (0,1) , 
and   are    arranged    so    that 


Rl    "    2^IR2    -   jl    *   ••■    ±    lR5    "    I'  (A"4) 


and   if   S    is    chosen   such   that   S    =    R.   with  probability    3/4    and 
S   =    R^   with  probability    1/4,    then    S   has   probability    density 
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2 
function    30x (1-x) (l-2x)        [3].      This    density    function,    it   was 

discovered,    requires    some    complex   programming   to   generate. 


E.       OTHER    DISTRIBUTIONS 

Let   Y  be    a   uniformly   distributed   random  variable    over   the 
interval    (0,1),    and    let 


y      =      F(x) 


/•y 


h(t)dt  (A-5) 


where  h(t)  =  1  for  0  t  1,  and  h(t)  =  0  otherwise.   Since 

x 

F(x)  px(x)  (A-6) 

> 

where  pv(x)  is  the  probability  density  function  for  X,  if 
x 

F(x)  or  p  (x) ,  the  desired  function,  are  known  explicitly, 
then  the  inverse  function 

x   =   F_1(y)  (A-7) 

will  generate  x  according  to  its  desired  density  function 
p  (x)  [3].   Most  often  pv(x)  will  be  a  polynomial  fitted  to 
some  statistical  data  which  represents  the  distribution  of 
the  variable  X.   Since  this  polynomial  is  normally  of  a  high 
order,  some  difficulty  can  result  from  this  method  in  that 
the  zeros  of  (A-6)  may  be  difficult  to  find. 
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APPENDIX    B 


COMPUTER    PROGRAMS 


C  COMPUTATION    OF    CUMULATIVE    0 1 ST^ I ^UTI CN    ^UNCTION    OF    3 

C  CB    ANGULAR     FPEOUENCY    OP    THE    SEP  US    R-L-C    CIRCUIT 

C  EXAMPLE.       PERFORMANCE    MEASURE     IS    AM    EXPLICIT    FUNCTION 

C  OP    THE    VARIABLES.        METHOD    IS    6    POINT    GAUSS    HPP^ITE 

C  QUADRATURE    SUM    APPROXIMATION    TO    EQUATION    (15)     TO 

C  OBTAIN    PROBABILITY    DENSITY    FUNCTION.       RECTANGULAR 

C  INTEGRATION     IS    THEN    USED    TO    OBTAIN    CUMULATIVE 

C  DISTRIBUTION    FUNCTION. 

COMMON  Z 

DIMENSION  t(6)fY(6),X(2) 

WRITE (6, 61 ) 

F  =  0. 
C      INPUT  GAUSS  HEPMITE  COEFFICIENTS  AND  ZEROS 

A(l)=.  453001E-2 

M2)  =  .  1570  6C7 

A(3)=. 7246296 

A(4)=A(3) 

A  (  5  )  =  A  (  2  ) 

A(6)=A(1) 

Y(1)=2.3506C5 

Y(2)=l. 335S49 

Y(3)=. 4360741 

Y(4)=-Y(3) 

Y(5)=-Y{ 2) 

Y(6)=-Y(l ) 
C      INPUT  WORST  CASE  FREQUENCIES  AND  COMPUTE  DELTA  Z  AND 
C      CISCRETF  FREQUENCIES. 

ZHI=.1879185E7 

ZLO=. 1405987c 7 

DELZ=(ZHI-ZL0)/19. 

Z  =  ZLO 
C  FOR    EACH    Z,    EVALUATE    PDF    USING    6    POINT    GAUSS    H5RMITE 

C  QUADRATURE. 

DC    3    IZ=1,2C 

P=0. 

DC    2    1=1.6 

DO    2    J  =  ]  ,6 

X( 1)  =  Y( J  ) 

X(2)=Y(  I  ) 

2  P=P+A( I  )*A(  J)*=UN(X) 

C  USE    RECTANGULAR     INTEGRATION    NOW    TO    OBTAIN    CUMULATIVE 

C  DISTRIBUTION. 

F=F+P*DELZ 
C  OUTPUT    CUMULATIVE    DISTRIBUTION 

WRITE (6 ,60) Z,F 
C  INCREMENT    Z 

3  Z=Z+DELZ 
STOP 

6C    cORMAT(  2E30.7) 

61     FCR"AT< » 1«. ////////, 15X. 'ANGULAR     FREQUENC Y • , 16X. 
l'DISTRI BUTICN"  ,//) 
END 
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c 
c 


c 

c 


c 
c 


FUNCTION  FUN(X) 

FUNCTION  SUBROUTINE  FOR 

QUADRATURE  SUM  ROUTINE. 

COMMON  Z 

DIMENSION  X(2),XO(3). 

FUN=0- 

CFFCK  FLAG  in    SEP.  tc 

IF( 11-121)1  ,3,  1 

PI=SORT( 3. 14159)**3 

P2=SORT(2.  ) 

INPUT  M«=£N  VALUES  Cc 

XC(  l}  =  1.5-3 

XO(2)=l .F-9 

X0(3)=l .53 

CALCULATE    STANDARD    DEVIATIONS 

CO    2    1=1.3 

SGMAU  )=XC(T)/30. 

SET    PL&G 

11=121 

CONTINUE 

PC    4    1  =  1  ,2 

"UMNORMALIZE"    VARIABLES    FOR 


USE    WITH    GAUSS    HERMITE 


SGMA(3),X1(3) 

CONSTANTS  HAVE 

COMPONENTS 


BEEN  CCVPUTEP 


OF  COMPONENTS 


AND  FUNCTION  INVERSE 
X1(I)=XU  )*P2*SGMAU)+X0U) 
X1(3)=Z*X1( 1)-1./(Z*X1(2) ) 
DRV=X1 (1)+1./(Z*Z*X1(2) ) 
COMPUTE  EXPONENT  Or  EXPONENTIAL 
F=(X1(3)-X0(2) )/ (P2*SGMA( 3) ) 

IP(F-13.6) 5,6,6 
EVALUATE  FUNCTION 
FUN=DRV*EXP(-F  )/<P2*°I*SGMA(3>  ) 

RETURN 
END 


COMPUTATION  CF  DERIVATIVE 


53 


ANGULAR     FREQUENCY 


DISTRIBUTION 


0. 1405 
0.14  3C 
0.1455 
0. 148C 
0.1505 
0.1530 
0.  1555 
0.15  80 
0.1605 
0.  1630 
1655 
1679 
1704 
1729 
0. 1754 
0.1779 
0.18C4 
0.1829 
0.18  54 
0.1879 


0 

0. 

0 

0 


C87E 
892E 
797E 
7C2E 
6C7E 
512E 
417E 
322E 
227E 
132E 
037E 
S42E 
847^ 
752F. 
£57H 
562E 
467E 
372F. 
277E 
182E 


C7 
07 

07 
07 
07 
07 
C7 
0  7 
07 
C7 
07 
07 
07 
07 
C7 
07 
07 
C7 
07 
07 


0, 

.9358540E- 

-06 

0. 

.2535363E- 

-04 

0. 

■  2094762E- 

-03 

0. 

2316327F- 

-02 

0. 

,  1303275E- 

-01 

0. 

,49607^°:- 

-01 

0. 

1364721E 

oc 

0. 

■2939899E 

00 

o. 

-4935315E 

00 

0. 

6891426E 

00 

0. 

.8430275E 

00 

0. 

.93148b9c 

00 

0. 

,97439031: 

00 

0. 

.9925191: 

00 

0. 

.  9980668E 

00 

0. 

,  999  48  28F 

0  0 

0. 

,9  998  8 13F 

00 

0, 

,  9999751E 

CO 

0. 

.9999889E 

CO 

o. 

.9999900E 

00 
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c 

cc 

c 

CB 

c 

EX. 

c 

VA 

c 

QL 

CC1 

CI 

WR 

c 

IN 

A( 

A( 

A( 

A( 

A( 

A( 

Y( 

Y( 

Y( 

Y( 

Y( 

Y( 

c 

IN 

c 

CI 

ZH 

ZL 

CE 

Z( 

CO 

IF 

Z( 

1  F( 

c 

FC 

DO 

OC 

DC 

X( 

X( 

X( 

Al 

2  F( 

c 

OU 

5  WR 

ST 

6C  FO 

61  FO 

1  «C 

EN 

MPUTATIO 
ANGULAR 
AMPLE  . 
RIA8LES. 
ADRATURE 
MMON  1(2 
MENSICN 
ITE(6,61 
PUT  GAUS 
1)=.1713 
2)=.36C7 
3)=. 4679 
4)=A(3> 
5)=A(2) 
6)=A(  1  ) 
1)=.9324 
2)=. 6612 
3)=. 2386 
4>=-Y(3) 
5)=-Y(2) 
6)=-Y(l) 
PUT  WQRS 
SCRETE  n 
I=. 18791 
0=. 14059 
LZ=(ZHI- 
1)=ZL0 

5  IZ=1, 

( IZ.cG.l 

IZ)=Z(  IZ 

IZ)=0. 

R  CACH  F 

2  1=1,6 

2  J=1.6 

2  K=l,6 

1)=Y(K) 

2)=Y( J) 
3)=Y( I ) 
=  A(I  )*M 

IZ)=F( IZ 

TPMT  CIS 

ITE(6,60 

IP 

P>4T< 2E3 

PMAT( ' 1« 

ISTRIBUT 

D 


CUMULATIVE  DISTRIBUTION  FUNCTION  OF 
OUENCY  Cc  THE  SERIES  R-L-CLC  CIRCUIT 
ORMANCE  MEASURE  IS  IWPLICIT  FUNCTION 
THCD  IS  6  °OINT  GAUSS  LEGENDRE 

APPROXIMATION  TO  EQUATION  (18). 
Z 
,Y(6) ,X(3),F(20) 


N  OF 
CR  = 
PERFORMANCE  MEASURE  IS  IWPLICIT  FUNCTION  OF 

ME 

SUM 
0),  I 
A(6) 
) 

S  LEGENDRE  COEFFICIENTS  AND  ZEROS. 

245 

616 

139 


695 
094 
192 


T  CASE  FREQUENCIES  ANC  COMPUTE  CELTA  Z  AND 

REOUENCIES. 

8  5  "7 

87E7 

ZLO/19. 

20 

)  GO  TC  1 

-D+DELZ 

pECU£NCYt  CCVPUTE  6  POINT  GAUSS  LEGENDRE  SUM 


J)*A(K) 

}+Al*FUN(X) 
TRIBUTION 
)Z(IZ) ,F{  IZ) 

0.7) 

.////////. 15X, 'ANGULAR 

ICN«  ,//) 


FREQUENCY* .  16X 
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FUNCTION  FUN(X) 
C      FUNCTION  SUBROUTINE  cOR  USE  WITH  G^USS  LEGENOPE 
C      INTEGRATION  SCHEME. 

CCMMGN  Z(2C) ,1 Z 

DIMENSION  XC(3),SGMA(3),Y(3) 

01  MENS  I  ON  X(3) 

PI=S0RT(2.*3.14159)**3 
C      INPUT  MEAN  VALUES  OF  COMPONENTS 

XC(1)=1 .E3 

XC( 2>=l.E-3 

X0<3)=1  .5-9 
C      CALCULATE  CCMP0N5NT  DEVIATIONS 

DO  4  1  =  1,3 

SGMA(  I)=X0(  D/30. 
C  "UNNORVALI Lc"    VARIABLES 

A    Y(I)=3.*SGVMI  )*X(I  )  +  XO(I) 

P=Y(1)/(2.*Y(2  )  ) 
C  COMPUTE     VALUE    GF    FER-CRMANCE    MEASURE    FOR    THESE 

C  COMPONENT    VALUES. 

ZZ=P+SORT( P*P  +  1 ./(Y(2)-Y(3) ) ) 
C  CHECK    TO    SEE    Ic     PERFORMANCE    MEASURE    IS    LESS    THAN 

C  REFERENCE    VA!  Up  . 

IF(ZZ.LE.Z(IZ)>CC    T G    1 
C  IF    NOT,     THE    INTEGRAND    IS    ZERO 

FUN=0. 

GC  TO  3 
C      Ic  IT  IS  COMPUTE  VALUE  OF  INTEGRAND 

1  F=0. 

DC  2  1=1 ,3 

2  F=F+X( T )**2 
FUN=27.*EXP(-4*5*F)/PI 

3  RETURN 
END 
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ANGULAR    FREQUENCY 


OISTRIRUTICN 


0. 14C5 
0.143C 
0.1455 
0.148C 
0.1505 
0.153C 
0.1555 
158C 
1605 
1630 
1655 
0.1679 
0.1704 
1729 
1754 
1779 
0.18C4 
0.1829 
0.1854 
0.1879 


0 
0 
0 


C8  7E 

892E 
7S7E 
7C2E 
6075 
512E 
4  17F 
322c 
227E 
132E 
037E 
S42E 
847-" 
752E 
£57£ 
562E 
467r 
372f: 
211  c 
182r 


C7 
07 
07 
CI 
C7 
07 
07 
07 
07 
C7 
07 
C7 
07 
C7 
C7 
07 
07 
C7 
07 
C7 


0.0 

0.10826 

0.40259 

0.61146 

0.51075 

0.22 150 

0. 83347 

0. 19609 

0.39814 

0. 58928 

0.79122 

0.90197 

0.93836 

0.96852 

0.98413 

0.98672 

0.98740 

C.9874^ 

0.98744 

0.98744 


3  2Lr- 

-05 

77-- 

-04 

64E- 

-03 

gQj-'. 

-02 

09:^- 

-01 

S2~- 

-01 

30: 

00 

95  = 

00 

13E 

00 

16; 

u  J 

74E 

00 

03r 

00 

29r 

00 

29:: 

00 

63^ 

00 

5  7:: 

0  0 

3  85; 

CO 

58  f 

CO 

59E 

00 
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c 
c 
c 
c 
c 
c 


c 
c 
c 


c 
c 


1 


60 
61 


COMPUTA 
CB  ANGU 

EXAMPLE 

PERFORM 

VALUES 

GENERAT 

CIMENSI 

DIMENSI 

WRTT£(6 

INPUT    W 

ZLO=.14 

ZHI=.13 

DELZ=( Z 

Z( 1)=ZL 

INITIAL 

F(1)=0. 

no   l   i  = 

Z(I)=Z( 

F(I )=0. 
TNJPIJT  " 
X0(1)=1 
X0(2)=l 
XC<3)=1 
DO  2  1  = 
SC-f'A<  I  ) 
INPUT  K 
IX(1)=1 
IX (2) =5 
I X  (  3  )  =  1 
SAMPLE 
DC  4  1  = 
GENERAT 
DC  3  J  = 
CALL  GP 
COMPUTE 
PP=X(1) 
ZZ=PP+S 
FIND  DI 
MEASURE 
DISTP IP 
DC  4  J  = 
IF(ZZ.L 
CONTINU 
DIVIDE 
D I  STRIP 
DO  6  J  = 
F(J)=-( 
OUTPUT 
WRITEC6 
STCP 
FOR^AT( 
FGRMAT( 
1 'DISTRI 
END 


TION    OF    CUM'JLATTVE    DI  STP !  P  UTI  CN    PUNCTICN    CF    3 
LAR    FREQUENCY    Oc    THE    S3TIES    R-L-C    CIRCUIT 
METHOD    IS    CRUDE    mqntE    CARLO,    WHERE    THE 
ANCE    MEASURE     IS    SAMPLED    10CC    TIMPS.       COMPONENT 
ARE    OBTAINED    FROM    GAUSSIAN    RANDOM    NUMBER 
OPS. 

CN    !X(3) 

GN    Z(20) ,X(3 >,X0(3) ,SGM A (3 > , F< 20 > 
,61) 

ORST    CASE    FREQUENCIES    AND    COMPUTE    DELTA    Z 
C598E7 
79185^7 
HI-ZLO/19. 
G 
IZE    CISTRIBUTICN    TO    ZERO. 

2,  20 

1-1 )+DELZ 

EAN    VALUES    OF    VARIABLES    AND    COMPUTE    DEVIATION. 

.E3 

•  E-3 

.E-9 

1,3 

=  X0(  D/30. 

ERNALS  FOR  RANDOM  GENERATION 

34^5 

4621 

74327 

PERFORMANCE  MEASURE  1000  TI^ES 

1,1000 

E    CO"cONENT    VALUES    FROM    GAUSSIAN    DISTRIBUTIONS. 

1.3 

N(  TX( J)  , SGMA (J)  ,X0( J)  ,X ( J) ) 

PER C0R VANCE  MEASURE. 
/'(2.*X(2I) 

QRT(PP**2+1°.  /<X(2)*X(3)  )  ) 
SCRETE  VALUES  OF  Z  FOr  v^HIC^  PERFORMANCE 

IS  LESS  THAN  OR  EQUAL  TO,  AND  INCREMENT  THE 
UTION  COUNTERS. 
1,20 
E.Z(J))F(J)=F(J)+1. 

DISTRIBUTION  COUNTER  BY  1000  TO  OBTAIN 

UTION 

1,20 

J)*l .E-3 

DISTRIBUTION. 

,60)Z(J ).F(J) 

2C3C7) 

• 1» ,////////,  15X, 'ANGULAR    FREQUENCY*  , 16X. 

BUT  I  CM'  ,//) 
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SUBROUTINE 
C  SUBROUTINE 
C      VARIABLES. 

A=0. 

DC  1  1=1,12 

CALL  UPM(IX,IY,Y) 

IX=IY 
1  A=A+Y 

V=(A-6. >*S+AM 

RETURN 

END 


GRNCIXt S,AMtV) 

TO  GENERATE  GAUSSIAN  DISTRIBUTED  RANDOM 


SUBPn'JTTNF 
C  SUBROUTINE 
C      VARIABLES. 

IY=IX*65539 

IF(IY) 1,2,2 

1  IY=IY+2147483647+1 

2  YPL=IY 
YFL=YFL*.46  56613E-9 
RETURN 

END 


URNU  X,I  Y,Y  =  L) 

TO  GENERATE  UNIFORM  DISTRIBUTED  RANDOM 
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ANGULAR  FREQUENCY 


DISTP IBUTION 


14C5 

1430 
1455 
1480 

15C5 
1530 
1555 
1580 
1605 
1630 
0.1655 
0.1679 
17C4 
1729 
1754 
1779 
18C4 
1829 
0. 1854 
0.1879 


0 
0 
0 
0 
0 
0 
0 
0 
0 
0 


930 
885 
79C 
695 
60C 
5C5 
41C 
215 
220 
125 
C30 
935 
840 
145 
650 
555 
460 
365 
27C 
175 


E  07 
-  07 
E  C7 
E  07 
E  07 

r  07 

S  07 

E  07 
E  07 
E  07 
F  C7 
E  07 
E  C7 
c  C7 
E  07 
r:  C7 
E  07 
E  07 
E  C7 
z    07 


0.0 

0.0 

0.0 

0.0 

0.6999999 

C.260000C 

0.7799995 

0.2120000 

0.3849999 

0.5929999 

0.75499*39 

0.8839999 

0.9489999 

0.9789^99 

0.9949999 

O.9090909 

0.9999999 

0.999O999 

0.9999999 

0.9999999 


r-02 
=  -01 
E-01 
c  00 
E  00 
E  00 
E  00 
E  CO 
F  00 
E  00 
E  00 
E  00 
E  00 
E  00 
00 
00 


E 
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c 
c 
c 
c 
c 


c 
c 


c 
c 


c 

CCMPUT^     "CCRNER    VALUES"    OF    I  NT E 

c 

^OUR    VALUES. 

DC    2    1=1,3 

X(I )=0. 

2    XKI  )=1.-X(  I  ) 

F1=FUN<  X  )+FU'](Xl  ) 

c 
c 
c 


ICO 
2C0 
3C0 

ACO 

500 
600 

7C0 

800 


C 

c 


COMPUTATION 
CB  ANGULAR 

PERFORMANCE 

VARIABLFS. 
CUBIC  ZERO 
CC^MCN  Z ( 20 
DIMENSION  X 
kRITE(6T61) 
INPUT  WORST 
DISCRETE  VA 
ZHI=. 187S18 
ZLO=. 140598 
DELZ=(ZHI-Z 
Z(  1)=ZL0 
CC  1  1=2,20 
Z(I)=Z( 1-1) 
DO  FOLLOWIM 
DO  999  IZ=1 
IN^UT  KERNA 
IX( 1)=49431 
IX(2)=39<+16 
IX(3)=17553 
SUM=0. 
CCM 
POU 
DC 
X(I 
Xl( 
Fl  = 

X(3)=l. 
XK31-0. 
F2  =  FUN( X  )  +  F 
X(  1)  =  1. 
Xl(  1)  =  0. 
F3=FUN(X)+F 
X(3)=0. 
Xl(3)=l  . 
F4=FUN(X)+C 
COMPUTc  I  NT 
AVG=(C1+F2+ 
NOW  GET  THE 
DO  7  JX=1T  2 
GENERATE  3 
DO  800  J=l, 
DO  6  1=1,3 
CALL  URN (IX 
I X  (  J  )  =  T  X  1  (  J 
IF(Y(1 )-Y(2 
IF(Y( 1)-Y( 3 
IF(Y(2 )-Y( 3 
IF(Y(1 )-Y(3 
X( J)=Y( 2) 
GO  TO  800 
X( J)=Y(1) 
GO  TO  8 CO 
X( J)=Y(3) 
GC  TO  SCO 
IF(Y(2)-Y( 3 
XI ( J)=l .-X( 
COMPUTE  REF 
SUM1=0. 
SUM1=SUM1+F 
SUM1=SUM1-F 
1X(2)*X(3J  ) 

SLM1  =  SUVJ,1-F 
SUM  1= SUM  1-^ 

SUM1  =  SU'U-P 

SLM=SUM+SUM 

COMPUTE  THE 

VALUE  Cc  Z 

P( IZ )=.C5*S 

F(  IZ)=F(  IZ  ) 


HF  CUMUL4T 

cREniJENCY  0 
MEASURE  IS 
METHOD  IS 

VARIANCE  ES 

),IZ 

(3),X1(3),Y 

CASE  FRcOU 
LUES  OF  Z. 
5^7 
7F7 
L0)/19. 


+  DELZ 

G  FOR  EACH 

,20 

LS  FOR  RAND 


IVE  DISTRIBUTION  FUNCTION  OF  3 
<=  SERIES  C-L-C  CIRCUIT  EXAMPLE. 

AN  IMPLICIT  FUNCTICN  OF  ITS 
TRIPLE  INTEGRAL  EXTENSION  TO 
TIMATOR  APPLIED  TO  EQUATION  (18) 

( 3), IX(3) ,1X1(3) ,F( 2C) 

ENCIES  AND  COMPUTE  DELTA  Z  AND 


REFERENCE 
CM  NUMBER 


VALUE    OF    Z 

GENERATION. 

GRAND    AND    GROUP    INTO 


UN (XI) 


LN(X1) 


UN (XI  ) 

TTAL    EST  I  MA 

=  3  +  F4)/8  . 

AVERAGE    OF 
C 
RANDOM    V4.RI 


TE 

20  PEFIN 
ABLES  WIT 


EVENTS 

H  6X( 1-X)  DISTRIBUTION. 


( J)  ,1X1  ( J)  ,Y(  I  )  ) 


} 

) )100, 10C0, 
)) 700, 1000, 
) )3C0, 10CC, 
) )500,1000, 


2C0 
500 
400 
600 


)  )4CC, 1CCC, 
J) 

I  NE  ME  NT 

UN(X )+FUN(X 
1*( l.-X(  1)- 

2*(X(3  )-X(l 
3*(  X(  2)-X(  1 
4*(X( 1)-X( 1 
1/<432.*X(I 

VALUE  DP  C 

U^+A VG 
/1.125 


6C0 


1) 
X(2)-X(3)+X(1)*X(2)+X(1)*X(3)+ 


)*X(3)-X( 

)*X( 2)-X( 
>*X(2)-X< 
>*X(2)*X( 

UMUL^TIVE 


2 )*X (3  )  +  X ( 1  )*X(2) ) 
3)*X(2)+X( 1 )*X(3) ) 
1)*X( 3)+X( 2)*X(3) ) 
3)*X1 (  1)*X1 (2)*X1  ( 

DISTRIBUTION  FOR 


3)  ) 

THIS 


61 


C      OUTPUT  DISTRIBUTION 
999  URITE(6,60)Z(IZ)  ,F(IZ) 
1000  STQ° 

60  FCRN"AT(2E30.7) 

61  FORMAT( 'l', ////////,  15X,«ANGULAF  FREQUENCY' ,16X , 
1»CISTRIBUTICN«.//1 

END 

FUNCTION  FUN(X) 
C      FUNCTION  SUBROUTINE  FCR  MONTE  CARLO  METHOD 

COMMON  Z( 20)  ,1  Z 

DIMENSION'  X0(3).SGMA(3) ,Y( 3). X(  3) 
C      CHECK  ^LAG  TO  DETERMINE  Ic  CONSTANTS  HAVE  BEEN 
C      FVALUATED. 

IF(II-114)1,3, 1 

1  PI=S0RT(2.*3.14159>**3 

C      INPUT  MEAN  VALUES  OF  COMPONENTS 

X0(1)=1  .E3 

XC(2)=l.r--3 

X0<3)=1  .E-9 
C  COMPUTE    DEVIATIONS     CF    COMPONENTS. 

DC    2     1  =  1  ,3 

2  SGMA(I)=X0(I)/30. 
C  SET    FLAG 

11=114 

3  G  =  0. 

C  B^-AK     INTEGRAND    INTO    8    PARTS    AND    SUM!    INDIVIDUAL 

C  FVALUATICNS. 

DO    6    1=1,2 
C  »UNN0RMALIZC"    VAQIA9LES     FOR    EACH    PART. 

Y ( 1 )  =  ( - 1 ) ** I *3 . *  SGM  * ( 1  )  *X ( 1 ) +X0 ( 1 ) 

DO    6    J=l,2 

Y(2)  =  (-l  )**J*3 .*SGVM2)*X(2)+X0(2) 

DO    6    K=l,2 

Y(3)=(-l )**K*3.*SGMA(3J*X( 3)+X0( 3) 

P  =  Y(1>/(2.*Y<2  )  ) 
C  COMPUTE    VALUE    Oc    PERFORMANCE    MEASURE 

ZZ=P+SQFT(P*P+1 ./(Y(2)*Y<3)  )  ) 
C  CHECK    VALUE    OF    PERFORMANCE    MEASURE    AGAINST    REFERENCE. 

IF(ZZ-Z(  !Z)  )4,4,6 
C  IF    LESS    THAN    OR    EQUAL    TO    REFERENCE,     EVALUATE    FUNCTION 

4  F  =  0. 

DO    5   l.=  l,3 

5  F=c+x(L)**2 
G=G+EXD(-4.5-F) 

C      IF  GREATER  THAN  REFERENCE.  SET  FUNCTION  ECUAL  TO  ZERO, 

6  CCNTINUE 
FUN=27.*G/PI 
RETURN 

END 


SUBROUTINE  «JRN(  I  X,  I  Y»  YFL  ) 

SUBROUTINE  CCR  GENERATING  UNITERM  RANDOM  NUMBERS. 

IY=IX*65539 

IF(  IYH.2.2 

IY=I Y+2147483647+1 

YCL=IY 

YFL=YFL*.4656613E-^ 

RETURN 

END 
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ANGULAR  CP?QU":MCY 


OTSTc  Ic'HTTOK< 


3. 

1405987" 

07 

0. 

143C892' 

07 

0. 

1455797" 

07 

0, 

14«0702" 

07 

0. 

15056C7? 

37 

0, 

1530  512= 

"\"7 

0. 

155541 7r 

07 

0. 

1580322* 

07 

0. 

1605227" 

07 

0. 

16  30132 1 

07 

0. 

16  5503  7- 

Z7 

0. 

]  6  79^42': 

0  7 

0. 

1 704  847" 

07 

0. 

172c>7^?,: 

r  7 

0. 

1754  65  7" 

^7 

0. 

1779562'-" 

07 

0. 

1804467  : 

07 

0. 

18  29372" 

07 

0. 

1854277F 

2  7 

0. 

1879182= 

07 

0.89 
0.89 
0.  39 
0.45 
0.26 
0.29 
0.11 
0.  17 
0.25 
0.30 
0.87 
0.92 
0.10 
0.  10 
0.  in 
0.10 
0.10 
0.1C 
0.  1C 
0.  10 


84364' 

84364- 

00373" 

55574' 

44415' 

9^560: 

55540: 

23926 

06310' 

67204' 
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